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1— I . 

^ ' Abstract. We report ab initio calculations of the melting curve and Hugoniot of molybdenum 

for the pressure range — 400 GPa, using density functional theory (DFT) in the projector 
augmented wave (PAW) implementation. We use the "reference coexistence" technique to 
overcome uncertainties inherent in earlier DFT calculations of the melting curve of Mo. Our 
calculated melting curve agrees well with experiment at ambient pressure and is consistent with 
shock data at high pressure, but does not agree with the high pressure melting curve from 
static compression experiments. Our calculated P{V) and T{P) Hugoniot relations agree well 
r~* . , with shock measurements. We use calculations of phonon dispersion relations as a function of 

pressure to eliminate some possible interpretations of the solid-solid phase transition observed 
in shock experiments on Mo. 



o 

^ I 1. Introduction 

' The melting curves of transition metals at pressures up to the megabar region are highly 

controversial, particularly for b.c.c. metals. Diamond anvil cell (DAC) measurements find that 
the melting temperature Tm increases by only a few hundred K over the range 1 - 100 GPa [HE], 
while shock experiments indicate an increase of several thousand K over this range [3l[ll[5]. There 
have been several ab initio calculations on transition- metal Tm(P) curves, and the predictions 
' agree more closely with the shock data than with the DAC data [U [71 [U [9]. A challenging 

. case is molybdenum, where there are very large differences between DAC and shock data [l2j . 

and where the shock data reveal two transitions, the one at high pressure (~ 380 GPa) being 
attributed to melting, and the one at low pressure (~ 210 GPa) to a transition from b.c.c. to an 
unidentified structure [3]. We report here on new ab initio calculations of T^{P) for Mo, and 
' on the P{V) and T{P) relations on the Hugoniot. We also report preliminary information that 

may help in searching for the unidentified high-P solid phase of Mo. 

We use density functional theory (DFT), which gives very accurate predictions for many 
properties of transition metals, including Hugoniot curves [Tn|. DFT molecular dynamics (m.d.) 
was first used to study solid- liquid equilibrium 12 years ago [llj, and several different techniques 
are now available for using it to calculate melting curves. In such calculations, no empirical 
model is used to describe the interactions between the atoms, but instead the full electronic 
structure, and hence the total energy and the forces on the atoms, is recalculated at each time 
step. There have been earlier DFT calculations on the melting of Mo, but the techniques used 
were prone to superheating effects [3 [9]. In the present work, we use the so-called "reference 
coexistence " techniques \19\ [T7\ [T^ , which does not suffer from this problem. 
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Figure 1. Comparison between PAW and FP-LAPW results for the GGA(PBE) and LDA(CA) 
approximations for E^c- Solid and dashed curves show GGA(PBE) and LDA(CA) FP-LAPW 
results, respectively; short-dashed and dotted curves show GGA(PBE) and LDA(CA) PAW 
calculations, respectively. Solid dots show experimental data [5j. 

Our work has several aims. First, we want to improve on the reliability and accuracy of 
the predicted melting curve of Mo obtained from DFT; second, we use DFT to predict the 
P{V) and T{P) relations on the shock Hugoniot; third, we want to identify the unknown solid 
phase of Mo observed in shock experiments. Our tests on the accuracy of DFT for Mo, and 
our extensive calculations of the Mo melting curve will be reported in detail elsewhere [TH], so 
we present only a summary here. However, our very recent calculations on the shock Hugoniot 
will be presented in more detail. These are important, because temperature is very difficult to 
measure in shock experiments [35] and DFT gives a way of supplying what is missing in the 
shock data. Our search for the unidentified solid structure of Mo is at the exploratory stage, 
but we present results for phonon frequencies as a function of pressure, which allow us to rule 
out some possibilities. 

In the following, we summarise our tests on the accuracy of the DFT techniques (Sec. [2]), 
and outline the reference coexistence technique. In Sec. [3] we present our results for the DFT 
melting curve and Hugoniot of Mo up to 400 GPa, and the study of the phonon frequencies. 
Discussion and conclusions are in Sec. HI 

2. Techniques and tests 

Our DFT calculations are performed mainly with the projector augmented wave (PAW) 
implementation [20], using the VASP code [2H I25j . since PAW is known to be very accurate. 
The main uncontrollable approximation in DFT is the form adopted for the exchange-correlation 
functional E^c- To test the accuracy of PAW, and the effect of E^c, we have compared our 
predictions for the P{V) relation of the b.c.c. Mo crystal against low-T experimental results 
(Fig. [1]). The pressure predictions from PAW using the Perdew-Burke-Ernzerhof (PBE) [28] and 
local-density approximation (LDA) [27J forms of E^c deviate by ~ 1.5 % in opposite directions 
from the experimental data, but we adopt the PBE form, because the deviations in this case 
are rather constant. The PBE results of Fig. [T]were obtained with 4s states and below in the 
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Figure 2. Comparison of calculated (curves) and experimental (dots) phonon dispersion 
relations of Mo at zero pressure. Experimental data are from Ref. [26j . 

core and all other states in the valence set. Inclusion of 4s states in the valence set makes 
no appreciable difference to the PBE results. We also tested the PAW implementation itself by 
repeating the P{V) calculations with the even more accurate full-potential linearized augmented 
plane- wave (FP-LAPW) technique [22l US |2l] , using the WIEN2k code [33j. As shown in Fig.lH 
PAW and FP-LAPW results are almost identical. Further confirmation for the accuracy of the 
PAW implementation and the PBE functional comes from our comparisons of the calculated 
phonon dispersion relations for the ambient-pressure b.c.c. crystal with experiment (Fig. [2j). 

The "reference coexistence" technique for calculating ab initio melting curves has been 
described elsewhere [18j . but we recall the main steps. First, an empirical reference model 
is fitted to DFT m.d. simulations of the solid and liquid at thermodynamic states close to 
the expected melting curve. Then, the reference model is used to perform simulations on large 
systems in which solid and liquid coexist, to obtain points on the melting curve of the model. 
In crucial third stage, differences between the reference and DFT total energy functions are 
used to correct the melting properties of the reference model, to obtain the ab initio melting 
curve. In the present work, the total energy function of the reference model is represented by the 
embedded-atom model (EAM) [SOlEl], consisting of a repulsive inverse-power pair potential, 
and an embedding term describing the d-band bonding. The detailed procedure for fitting Urei 
to DFT m.d. data will be reported elsewhere [18|, but we note that we needed to re-fit the 
model in different pressure ranges. 

The simulations on solid and liquid Mo in stable thermodynamic coexistence using the fitted 
reference model employed systems of 6750 atoms, and we checked that the results do not change 
if even larger systems are used. The protocols for preparing these simulated systems, and for 
achieving and monitoring stable coexistence were similar to those used in earlier work on the 
melting curve of Cu [TTj. The procedures for correcting the melting curve of the reference model, 
which depend on calculations of the the free energy differences between the reference and DFT 
systems, have been described and validated in earlier work [TH]. 
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Figure 3. Calculated ah initio melting curve (filled circles and solid line) of this work 
compared with previous results: generalized pseudopotential calculations of Moriarty [6] (dotted 
line), dislocation-mediated models of Belonoshko et al. [9] (long-dashed line) and Verma et 
al. [29] (dashed-dotted line); experimental shock- wave [3] and DAC [Ij measurements are shown 
with empty squares and triangles, respectively. Filled and inverted-empty triangles show solid 
and liquid ah initio molecular dynamics calculations of Belonoshko et al. [5] , respectively. Empty 
circles show results of this work obtained with the EAM model without free-energy corrections. 



3. Results 

3.1. Melting properties 

Fig. [3] shows the melting curve of our reference EAM model, and the melting curve obtained 
from this by correcting for the differences between the DFT and reference total-energy functions; 
earlier ah initio-hased calculations of the Mo melting curve due to Moriarty and to Belonoshko 
et al. are also indicated [9]. We also show points on the melting curve from DAC and shock 
measurements [H H] . The differences between reference and corrected melting curves are only 
a few hundred K, so that the corrected curve should be very close to the melting curve that 
would be obtained from the (PBE) exchange-correlation functional if no statistical-mechanical 
approximations were made. Because we avoid approximations of earlier ah initio-hased work, 
our present results should be a more accurate representation of the DFT melting curve. Up to 
100 GPa, the differences between our results and earlier DFT work are rather small, and we 
confirm that DFT gives a much higher melting slope than that given by DAC experiments. We 
obtain an accurate value of dT^/dP at P = by fitting our melting curve to the Simon equation 
Tm = a(l + P/b)", with o = 2894 K, 6 = 37.2 GPa, c = 0.433. The resulting P = value of 
Tin = 2894 K is close to the accepted experimental value Tm = 2883 K. Our dT^/dP value of 
33.7 K GPa~^ at P = agrees with an older experimental value of 33.3 K GPa~^ [32]. Our DFT 
melting curve is consistent with the point obtained at P ~ 370 GPa from shock measurements. 
However, we stress that the temperature of the "experimental" point was not measured, but 
estimated by considering superheating corrections to the shock- wave data [121113]. Because of 
this, it is important to try and corroborate the estimated experimental temperature, and this 
can be done by ah initio calculations, as we explain next. 
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Figure 4. Left: Pniyn) ah initio relation on the Hugoniot of Mo up to pressure of 400 GPa; 
dots reproduce experimental data of Refs. Right: T(Pi{) ah initio relation on the Hugoniot 

of Mo up to pressure of 400 GPa; dots reproduce experimental data of Refs. [H [5] and triangles 
the same results but corrected for superheating effects as given by Ref. [E]. Error bars show 
the uncertainties in pressure and temperature. 



3.2. Hugoniot curves 

The pressure Phj volume Fh and internal energy Ey{ in the shocked state are related to the 
initial volume Vq and internal energy Eq by the well-known Rankine-Hugoniot formula [34] : 

^Ph {Vq - Vn) = En-Eo. (1) 

Since the internal energy and pressure are given in terms of the Helmholtz free energy 
F hy E = F — T{dF/dT)v and P = —{dF/dV)T, we can calculate the Hugoniot from 
our DFT simulations, provided we can calculate F as a function of V and T. So far, 
we have done this only for the b.c.c. crystal in the harmonic approximation, in which 
F{V,T) = Fperf(V^,T) + -Fharm(^i Here, -Fperf is the free energy of the rigid perfect crystal, 
including thermal electronic excitations, and -Fharm is calculated from the phonon frequencies 
ujqs (q the wavevector, s the branch). We calculate i*harm in the classical limit, in which 
-^harm = SksT ln{ptiu}) per atom, with (3 = l/k^T and Co is the geometric average of phonon 
frequencies over the Brillouin zone. The methods used for to calculate Fperf (V^,T) and the 
frequencies ujqs were similar to those used in our earlier work on Fe (Ref. [IS]). For a set of 
temperatures, we calculated -Fperf at a set of volumes, and fitted the volume dependence with 
a third-order Birch-Murnaghan equation [36]. The temperature dependence of the coefficients 
in this equation were then fitted with a third-order polynomial. The phonon frequencies were 
calculated at 12 volumes in the range 15.6 — 9.2 A^/atom, as explained elsewhere [18j. The 
volume dependence of the average ui was then fitted with a third-order polynomial. 

To obtain Ph(^h) and T{Pii) from our fitted free energy, for each value of Vu we seek the T 
at which the Rankine-Hugoniot equation is satisfied, and from this we obtain Ph- For Vq and 
£^0) we used values from our GGA calculations; we checked that use of the experimental Vq made 
no significant difference. Comparison of our calculated PuiVn) with measurements of Hixson 
et al. [U [5] (left panel of Fig. U]) shows excellent agreement. In the right panel, we compare 
our T(Ph) with both uncorrected results of Hixson et al. and also with results corrected for 
superheating, and our results confirm their temperature estimates. 
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Figure 5. Comparison of the calculated phonon dispersion relations of Mo at different pressures: 
GPa ( ), 77 GPa ( ), 136 GPa ( ), 274 GPa ( ) and 400 GPa ( ). 

3.3. Solid-solid phase transition 

Efforts have been made to identify the high-P/high-T structure of Mo indicated by shock 
experiments [4J. Hixson et al. used their theoretical prediction of b.c.c. ^ h.c.p. phase transition 
at P ~ 320 GPa and T = K to suggest the h.c.p structure. However, later calculations locate 
this transition at higher pressures (420 < P < 620 GPa), and temperature stabilisation of h.c.p. 
over b.c.c. below melting seems improbable. Furthermore, there are recent claims that under 
pressure b.c.c Mo transforms first to f.c.c. rather than h.c.p. [37]. Very recently, Belonoshko 
et al. [9J reported ah initio simulations on the f.c.c and A15 structures at temperatures and 
pressures near those of interest (namely, P ~ 210 GPa and T ~ 4000 K), and concluded that 
both structures are unlikely high-T phases of Mo. We have performed ah initio calculations 
similar to those of Belonoshko et al. on the h.c.p. and u} structures [P ^ 250 GPa and 
T = 5000 K). We find that temperature neither favours any of them respect to b.c.c, thus they 
must be excluded as well. Recently, a 2nd-order phase transition from cubic to rhombohedral 
has been observed in vanadium at P = 69 GPa [38]. This appears to be related to an earlier 
ah initio prediction of a phonon softening in V [39j. This suggests that a similar structural 
transition might occur in Mo. We have used our calculated phonon dispersion relations of Mo 
over the range 0-400 GPa to test this. Fig. [5] depicts our results at 0, 77, 136, 274 and 400 GPa, 
but we see that no phonon anomaly like that reported for V is observed at any wavevector. This 
indicates that we should rule out structures based on small distortion of b.c.c. 

4. Discussion 

An important outcome of the present work is improved DFT calculations of the melting curve 
of Mo over the pressure range — 400 GPa. In particular, our techniques allow us to avoid the 
superheating errors which appear to affect an independent recent DFT study of Mo melting [9]. 
The accuracy of our calculations is confirmed by the very close agreement with experiment for 
the melting temperature and the melting slope dT^/dP at ambient pressure. The results 
fully confirm that the increase of Tm by ~ 2000 K over the range — 100 GPa predicted by DFT 
is about 10 greater than that deduced from DAG measurements. A second important outcome is 



that our calculations of the temperature along the shock Hugoniot support earlier temperature 
estimates based on experimental data but corrected for possible superheating effects [12]. This 
allows us to compare more confidently the point on the melting curve at P ~ 380 GPa derived 
from shock data with our predicted melting curve, and we confirm that at this pressure is 
ca. 8650 K. This is far above any reasonable extrapolation of the DAC data. Concerning the 
search for the unknown crystal structure of Mo indicated by shock experiments to exist above 
ca. 220 GPa, we have been able so far only to rule out some possibilities. Our calculations of 
the phonon dispersion realations in the b.c.c. structure over the range — 400 GPa reveal no 
softening of any phonons, and no indication of any elastic instability. This means that the new 
crystal structure does not arise from small distortion of b.c.c. This is interesting in the light of 
the recent discovery of the elastic instability of b.c.c. V above P ~ 70 GPa, predicted initially 
by DFT, and observed very recently in x-ray diffraction experiments [38j. It seems that the 
structural transition in Mo is of a different kind. 

The large conflict between the melting curve of Mo derived from DAC measurements on 
one side and from shock experiments and DFT calculations on the other side must be due 
either to a misinterpretation of the DAC data or to a combination of serious DFT errors 
and misinterpretation of shock data. Given the accuracy of DFT that we have been able to 
demonstrate (low-temperature P{V) curve, Hugoniot P{V) curve, ambient-P phonons), we 
think there is little evidence for significant errors in DFT, which is also in good accord with 
shock data. A possible explanation might be that the large temperature gradients and non- 
hydrostatic stress in DAC experiments might give rise to flow of material giving the appearance of 
melting, even well below the thermodynamic melting temperature. We also note recent evidence 
that temperature measurement in DAC experiments may be subject to previously unsuspected 
errors [lOj, though probably not of the size needed to resolve the conflict by themselves. 
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